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Abstract 

We study the effect of a constant electrical field applied on vicinal surfaces such as the Si(lll) 
surface. An electrical field parallel to the steps induces a meandering instability with a nonzero 
phase shift. Using the Burton-Cabrera-Frank model, we extend the linear stability analysis per- 
formed by Liu, Weeks and Kandel (Phys. Rev. Lett. 81, p. 2743, 1998) to the nonlinear regime for 
which the meandering amplitude is large. We derive an amplitude equation for the step dynamics 
using a highly nonlinear expansion method. We investigate numerically two limiting regimes (small 
and large attachment lengths) which both reveal long-time coarsening dynamics. 

PACS numbers: 66.30.Qa, 47.20.Hw 
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I. INTRODUCTION 



Stepped crystal surfaces exhibit a number of different morphological instabilities likely to 
play an important role during crystal growthi^^iMi^. Furthermore, the ability to control the 
growth of faceted stepped crystal surfaces may be of considerable importance when manu- 
facturing electronic and optoelectronic devices^. These morphological instabilities occur not 
only during growth and evaporation but also under the influence of an electrical field, as on 
the well-studied Si(lll) surface of a semiconductor. Surface electromigration instabilities 
may also arise in metals, where they are an important source of failure in microelectronic 
devices at metal-metal interfaces and also an interesting tool for pattern formationi2*S. One 
of the most studied instability, known as step bunching, arises on the Si(lll) surfaces from 
the biased diffusion (drift) of adatoms under the influence of an external driving force such as 
an electrical constant fieldiSiiiJ^iiSiii. Step bunching is a one- dimensional instability which 
has been explained within the framework of the Burton-Cabrera Frank equationsA^ in terms 
of displacements of steps and terraces. Recent experimental and theoretical studies of step 
bunching revealed several difficulties, like the complex role of step transparency, the Ehrlich- 
Schwoebel barriers, the effect of substrate temperature, and the variations of the adatom 
mobility with the distance to the step3i^*i2*i2^. 

In the present study, a constant electrical field is applied along the mean step direction 
of a train of synchronized steps (all identical up to a constant phase-shift). An experimen- 
tal study of a comparable system was recently reported^^i and it was shown in this work 
that a two-dimensional step meandering instability takes place. The linear analysis of this 
problem was previously performed by Liu and co-authors who predicted the occurrence of 
synchronized meandering2i. We perform here a nonlinear analysis of this instability in order 
to describe the long time behavior of the in-phase meandering mode. In particular we show 
the appearance of the coarsening regime in which steps undulations increase. This paper is 
organized as follows. In the next section, we present a model based on the Burton-Cabrera- 
Frank equationsi^. In the third section, we perform the linear analysis, which serves as a 
basis for the nonlinear analysis. A general nonlinear evolution equation including the effects 
of the repulsive step-step interactions is derived in section IV. The results of numerical simu- 
lations of this nonlinear equation are presented and discussed in section V, while conclusions 
and perspectives are postponed to section VI. 
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Before presenting our model, we shortly review previous work concerning nonlinear equa- 
tions for the time evolution of synchronized steps. The step meandering instability was 
originally predicted theoretically by Bales and Zangwill^^ for a vicinal surface under growth. 
Its origin is the asymmetry between the descending and ascending currents of adatoms. As 
shown by Bales and Zangwill, a straight train of step may become morphologically unsta- 
ble during MBE growth if the kinetic attachment at the steps is asymmetrical: this is the 
Ehrlich-Schwoebel effect. It was shown that the most dangerous mode corresponds to a zero 
phase-shif1)2^. Nonlinear extensions of this work have shown that the meander evolution 
can be described by amplitude equations displaying diverse behaviors. Close to the in- 
stability threshold, starting from the Burton-Cabrera-Frank (BCF) model, it was proved^^ 
that the step position in the presence of desorption (evaporation) obeys the Kuramoto- 
Sivashinsky equation. The ultimate stage of this dynamics is thus spatio-temporal chaos. 
In the case of negligible desorption with strong or moderate Ehrlich-Schwoebel effect, it 
was found that the step amplitude obeys a highly nonlinear equation . ^^i^^i^^ This equation 
cannot be derived from a weakly nonlinear analysis but is based on the assumption that 
the slope of the steps is of order unity. Instead of spatio-temporal chaos, a regular pattern 
is revealed: the lateral modulation wavelength is fixed while the transverse amplitude of 
the step deformation (meandering amplitude) increases. Elasticity or diffusion anisotropy 
can also influence the meander dynamics.— It was recently shown in context of the step 
meandering instabilities during growth on Si(001)^2i21 that the nonlinear dynamics is driven 
by a conserved Kuramoto-Sivashinsky equation. This equation was already mentioned in 
Ref.-^ on the basis of symmetry arguments but was not derived there, because of a different 
scaling of the Ehrlich-Schwoebel effect. Step meandering was also studied in the context of 
electromigratioii22i2^iM using analytic linear analysis and kinetic Monte Carlo simulations. 
It is therefore of importance to extend the work of Liu et and to develop analytical tools 
to describe the nonlinear regime of the meandering instability. 
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FIG. 1: Schematic representation of a small portion of a vicinal surface showing three steps. On 
the top view (above), the electrical field E is represented. On the side view (below), two attachment 
mechanisms are illustrated. 

II. MODEL 

A. Validity range and notations 

The geometry of the problem is sketched in Fig. (^. Initially, all the step edges are 
directed along the X axis, and equidistant from a distance Lq. A constant electrical field E 
is applied in the positive X direction. To investigate the resulting meandering instability, we 
use a two-dimensional version of the BCF model. The terraces are numbered sequentially 
in the step-down direction. In our notations, the n-th terrace is bordered by the two steps 
numbered n and n + 1. To distinguish between quantities defined anywhere on the terrace 
and quantities defined at steps only, we will use an upper index n in the first case and a 
lower index n in the second. For instance, the adatom concentration on the n-th terrace is 
denoted C", while the equilibrium concentration at step n reads C^'^. 

In practice, it is usually assumed that the concentrations are not explicitly time- 
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dependent, so that C" = C"(X, F). This quasi-static approximation is vahd provided 
that the diffusion coefficient of the adatoms on the terraces is sufficiently large that 
diffusion takes place on time scales shorter than those for step motion. The diffusion bias 
introduced by the electrical field can be quantified by the ratio of the thermal energy fc^T 
to the electrical energy \Z*e\ElE- Balancing the two terms defines the electrical length as 

w 

where ks is the Boltzmann constant, T the absolute temperature, Z* the effective atomic 
charge number, and e the electron electrical charge. 

For the sake of simplicity, we set both deposition and evaporation of adatoms to zero here, 
whereas experiments are usually performed with a small but nonzero net flux. Introducing 
both effects in our model is straightforward and would not affect qualitatively the results 
obtained within the zero flux assumption. 

On the n-th terrace, the quasi-static biased diffusion equation reduces to 

Ds{dxx + 5ry)C" - {Ds/lE)dxC'^ = 0, (2) 

and the adatom flux is 

r = Ds{l/iE-dx,-dY)C\ (3) 

The boundary conditions for Eq. Q are obtained by writing mass conservation at all 
the points R„ = (X, 1^) and Rn+i = {X, l^+i) located on both edges of the terrace. In 
the present model, we assume that the adatom attachment/detachment kinetic coefficients 
are the same on the upper and lower side of a given step, = = v. We further re- 
strict ourselves to temperature ranges where direct mass exchange between adjacent terraces 
(transparency) can be neglected. At R„, the boundary condition thus reads 

J"(R„) ■ u„ = -z/[C"^(R„) - Q"], (4) 

where u„ represents the normal unit vector pointing in the step-down direction. Alterna- 
tively, at R„+i, the second boundary condition is 

J"(R„+i) ■ u„+i = z/[C"(R„+i) - (5) 

Writing mass conservation at any point R„ along step tt,, we obtain the normal (component 
along u„) step velocity, 

K = ^^.z/[C"(R„) + C"-i(R„,) - 2^"] . (6) 
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In this equation, fi^ is the adatom area and we neglect adatom diffusion along the step. 



B. Nondimensional version of the governing equations 

The terrace width Lq provides a natural length scale for the problem. Another possibity 

is the electrical length ^E defined in Eq. (Q). However, since is likely to diverge as the 

electrical field goes to zero, Lq is prefered. Setting 

X Y \ ij, ^ 

a; = ?/ = l-> - = c = — , (7) 

where CqQs is the fraction of adsorption sites occupied by adatoms, we get the nondimen- 
sional form of Eqs. (j2El)- The quasistatic diffusion equation reads 

{d.. + dyy - vd.)c^ = 0, (8) 

and the adatom flux is 

r = (r/-9.,-9,)c". (9) 
In this equation, the nondimensional flux vector is defined as 

j" = 7^ J", (10) 

SO that the physical time is rescaled by the characteristic time 

Using both the time and space scale factors, we obtain the first, 

r(rO-u„ = -p[c"(r„)-C], (12) 
and the second boundary condition, 

j"(r„+i) ■ u„+i = p[c"(r„+i) - cl%^]. (13) 
The nondimensional number, 

-f. (") 

inversely proportional to the characteristic length d = Dg/i', indicates which mechanism 
between diffusion and attachment governs the time evolution of the steps. The normal 
velocity of step n now takes the form, 

Vn = ap[c"(rO + c"-i(r„) - 20^"] , (15) 
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where 



^ = 72 (16) 



is the ratio of the two basic areas in the problem. 



C. Equilibrium concentration 

We now derive a detailed expression of the equilibrium concentration c^'^ at step n. A 
quite general form is 

C = C,tVCo = exp(^)=l + ^ + ... (17) 

We will use the thermal energy ksT as the energy scale, and define the nondimensional 
chemical potential as 

M , , 

Within the nondimensional description presented in previous section, the position of step n 
is represented by a function yn{x), and we define the relative position as 

Cn{x) = yn{x) - n. (19) 

Following Paulin and coworkers,^^ we introduce a nondimensional free energy functional for 
step ra, 

fn = + fl (20) 

The first term is due to the step stiffness. 



= ds. (21) 

J n 

Here ds is the integral of the curvilinear abscissa along the whole step n (total step length) 
and (3 = B{LQ/kBT), where B is the step-stiffness of the material. The second term sums 
up the step-step repulsive energies assumed to vary as the inverse square distance. 



ii^Y+i^y (22) 
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1 \ 2 / 1 \ 2 



n 



where a = A/lksTLQ), and A is the step interaction coefficient. As shown in Fig. Q, the 
lengths and l~ are the shortest distances between steps {n,n + 1) and steps {n,n — 1). 
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FIG. 2: Shortest distances between a given step and its two closest neighbors, in the general case. 
The tangents to steps n — 1 and n + l are drawn at two points having the same abscissa x. 

Thus the previous relation gives only an approximate value of the total repulsion energy. 
Finally, the chemical potential is obtained by a functional derivation of the free energy. 



III. LINEAR STABILITY ANALYSIS 

Repulsions between steps prevent them from intersecting one-another and edge stiffness 
limits their curvature. However, the possibility that the shapes of two consecutive steps are 
weakly correlated in phase or amplitude remains open. Since the general problem is quite 
difficult to solve in practice, we limit the present study to the simple case of a synchronized 
train of steps. 

Starting with straight steps, separated by a unit distance, we introduce an harmonic 
perturbation of amplitude e ^ 1, wave number q, and phase-shift 0, 




(23) 




(24) 
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Looking for solutions of the nondimensional diffusion equation under the form 



c^{x,y) = l+cl{y)Ux), (25) 

we obtain 

cl{y) = A^e'y + B'^e-'y, (26) 

with 



^/ q"^ + irjq. (27) 



To derive the dispersion relation, we first express the chemical potential using Eqs. 
In practice, the step curvature is small, and an accurate approximation of the step-step 
distances is given by 

/+= 1 + I- l + Cn-g ^ ^28) 

To the leading order in the perturbation amplitude e, we find the chemical potential from 
Eqs. T MMI^ . 

= (^Cn9iq,(l)), (29) 

where 

g{q,(f)) = {a + (3) q^ + 6a{l- cos (f)). (30) 
Introducing this result in Eqs. p2|lHpi5|l . the following dispersion relation is finally obtained, 

uj{q, 0) _ 9 J sin + a(cos - cosh r - ^ sinh r)g{q, 0) 
2ar ~ (l + ^)sinhr + 2^coshr 

This equation is similar to the dispersion relation derived in Ref.— . In the remaining of this 
section, we will neglect the step-step interactions, a = 0, to avoid unnecessarily complicated 
relations. Keeping the actual value of a would only introduce small quantitative changes. 
The growth rate is defined as r(g,0) = Re(co'). A meandering instability is thus expected 
for positive values of F. Experimentally, the electrical field is a weak perturbation, so that, 
according to Eq. (jZj), the parameter 7] takes small values (~ 10~^ — 10"'*). We thus expect 
the wave numbers of the corresponding instabilities to verify g <^ 1. As we do not know 
a priori the relative magnitude of rj and q, we first use a general expansion in which both 
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FIG. 3: Plot of the growth rate r(g, 0*) = Re{uj) obtained from the hnear stability analysis Eq. 
(|3H) as a function of the perturbation wave number q. 



quantities are small and considered equivalently. This leads to 



r 




2 2p 



2 



— COS 



0) 



+ 



(32) 



Although the two last terms are of higher order, we nevertheless keep them since the lowest 
order terms vanish when = 0. A first necessary condition for a maximally unstable mode, 
{dT/d(f))g = 0, gives the following phase-shift: 



The corresponding growth rate, r(g, 0*), is plotted as a function of the wave number q in 
Fig. ©: a meandering instability arises as soon as the electrical field is nonzero.— For q ^ r], 
it is possible to find relations between q and t] in different ranges of wave numbers. The 
most unstable mode q = qmax is obtained by introducing the second condition (dT/dq)^ = 0. 
Together with Eq. (jHHjl . we get 




(33) 




(34) 
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and the absolute maximum of the amphfication factor is 



/5p(2 + p) 



Finally, the marginal mode g = go is deduced from the condition r(g, 



*HiVf^Vfel. (36) 



The scalings of qmax and go with r] thus suggest that the range of unstable modes is quite 
large. The limit of weak electrical fields (r^ <^ 1) is relevant for the experimental work re- 
ported in Ref."2ii, in which step meandering is observed. In addition, the maximum growth 
rate being small, the situation is favorable for a nonlinear analysis of the meandering insta- 
bility which is presented in the next section. 

IV. NONLINEAR ANALYSIS 
A. Local coordinates 

As illustrated in Fig. (jlj), we consider the case of steps which are all identical up to 
a translation in an oblique direction z rotated by an angle Q with respect to axis y. The 
amplitudes of two successive steps have thus the following property: 

Cn+i(a; - tan6') = Cn(a;). (37) 

In the linear regime defined by Eq. (j23), we have 

~ tan-^((/)/g). (38) 

Since we want to explore nonlinear dynamics, the amplitude of the meanders may reach 
values of order unity for which Cra(2;) is no more a single- valued function in the original frame 
of reference (x, y). For this reason, our nonlinear model makes use of a non-orthogonal frame 
of reference (x,5), defined as 

X = X ^ ytanO, 

~z = (39) 
cosfe' 

as shown in Fig. (@]). With this change of coordinates, the step shape function becomes 

m = (40) 
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FIG. 4: A set of steps identical up to a translation along the z axis. The local variables used in 
the nonlinear analysis are z = y/{cos6) and x = x + ytanO. 

where the n index can be omitted because all the steps are identical in the new frame. We 
further define a local frame of reference {x, ip) by moving the x and z axes along the step, 

X = X, 

^ = i_^(5)=5_^(;^). (41) 

In the local frame, the partial derivatives transform as 

dy = tune + - tan^ d^^] d^, (42) 

Vcosy / 

where d^C, = d^/dx- The second derivatives are derived from these expressions. 
B. BCF equations for the local coordinates 

Introducing the relations for the partial derivatives into Eq. (jH)), one gets the quasi-static 
diffusion equation for the local coordinates (Xyi^): 

+ 2{sme - d^Odx^ - Vcos' e d^]c{x,^), (43) 

where 



p{x) = V(l - sin^ d^O' + {cose d^O'. (44) 
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Note that the step index n is purposely omitted because of translational invariance. It is 
easier to express the vectorial quantities in the base of the two unit vectors of the initial 
frame {x,y). The adatom flux reads 



j = (jxjy) = c{x, '^)y-d^ + {dx09^ + V , 

— tan^^ ~~n9^ + i^xO^^^^ (45) 

COS u / 

and the unit normal vector to the step, 

u = i - COS0 (9^^, 1 - sin^ d^^^ . (46) 
The two boundary conditions take on very simple forms, 

j ■ u = — p(c — c'^'^) at ip = 0, 

j-u = +p(c-c^«) at ^ = ^. (47) 

cosy 

The expression of the local curvature is needed to complete these boundary conditions. We 
obtain 

, , cos6'„ ^ 
<X) = -^d^x^- (48) 



The normal velocity is deduced from Eq. (|T^. 



V{X) = = CTp 

p 



c(x,0) + c(x, -^)-2c'%x) 
V cos / 



(49) 



cos 9^ 

where c^''(x) — 1 + fJ'ix)- An expression for the chemical potential is obtained through the 
functional derivative of the free energy functional, Eqs. ()20II23|1 . In the oblique frame of 
reference, we have now, 

and the curvilinear length element is ds = p dx- As illustrated in Fig. (0) the shortest 
step-step distances are defined in a slightly different way in this frame: the tangents to the 
two adjacent steps are drawn at a given value of x. Adapting Eq. pHj) to this new definition, 
we obtain 

/+ = -s^^il — i and /- = -i^^Li — ^, (51) 
P+ P- 

where 

P± = ^{l-sm9d^^±y + {cos9d^^i)'^. (52) 
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FIG. 5: Shortest distances between a given step and its two closest neighbors, in the local frame. 
The case of translational invariance along z is represented here. The tangents to the adjacent steps 
are drawn at two points having the same abscissa x. 



Note that it is necessary to keep the amplitudes ^, and ^+ of three successive steps to 
perform the functional derivation. After derivation, we set ^ = so that = l_ = I, and. 



a 



+ (^x^ - sin I 



(53) 



P cos^ 6 

The chemical potential sums up the contributions of the step stiffness and of the step-step 
interactions. 



C. Small parameter expansion 

The aim of this paragraph is to establish a nonlinear equation for the time evolution of 
a step. 



1. Scaled variables 

In the linear analysis presented above, we have shown that the wave number go ~ V^^'^ is 
small as compared to unity in the limit of a weak electrical field. We thus introduce a small 
parameter e <c; 1, such as 

V = (54) 
and rj2 is of order unity. As a consequence, we define the slow space variable 

X = ex- (55) 

Note that the slow x variable used hereafter differs from the fast x variable discussed in 
section II. This should not introduce any confusion, since only the new x appears in the 
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following. At the marginal wave number q = qq, the linear analysis results of Eqs. (jS^l 
Sni IHH|) give the following relation between the inclination angle 6 and the non-dimensional 
number p, 

^^^tan-iff^)^Y (56) 



P 

Since p can be large or small depending on the parameters, 6 can take arbitrary values. The 
boundary conditions given in Eq. (jTfjl are applied ai ip = and ip = 1/ cos6, so that the 
space variable ip is simply equal to z. Accordingly, one defines the meander amplitude as 
h{x) = ^(x) and the normal velocity as v{x) = v{x)- 

2. Order by order expansion 

The unknown variables, concentration and shape function are expressed as power expan- 
sions of the scaling parameter e, 

c{x, z) = co(x, z) + eci(x, z) + e'^C2{x, z) + . . . 
h{x) = e'^h_i{x) + e°ho{x) + ehi{x) + . . . 

v{x) = e^V3{x) (57) 
Introducing this expression of h{x) in Eq. (jH)), the following development is found for p, 

p{x) = po{x) + epi{x) + e^p2{x) + . . . , (58) 

with 



Po{x) = - sine d^h^if + (cos^ d^h^if. (59) 
We obtain in a similar way the equilibrium concentration 

c^'i^x) = 1 + ecf (x) + e^cl\x) + . . . , (60) 

where 



Ci (x) = —crcosO 5 — 

Po 



(3 + apl(2pl - cos^ e'^y (61) 



We now solve order by order the nondimensional equations obtained by introducing the 
scaled variables defined above into Eqs. ()43l 1^ 1471 . The results obtained at order i are 
used to derive the equations at order i + 1. 

order 
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The diffusion equation reduces to 

pI dzzCo = 0. 

We look for solutions of the form 

Co{x,z) = ao{x)z + bo{x), 

which imposes 

pI Oo - ppo cos 9 [bo - l] = 0, 

Po «o + PPo COS 6 [bo - l] + ppo do 

for the boundary conditions, so that, 

ao{x) = 0, 

bo{x) = 1, 

Co{x,z) = 1. 

At this order, the velocity is found to be zero. 

order 1 



Diffusion equation: 



Solution: 



Boundary conditions: 



Solution: 



pI dzzCi = 0. 
Ci{x, z) — ai{x)z + bi{x). 

pI oi - PPo cos 6 [bi - ] = 0, 
Po «i + PPo cos 9 \bi - ] + PPo ai 

ai{x) = 0, 
h^[x) = cl%x), 
ci{x,z) = cf(a;). 



At this order, the velocity is found to be zero. 
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order 2 

Diffusion equation: 

PO dzzC2 = 0. 

Solution: 

C2(x, z) = a2{x)z + b2{x). 

Boundary conditions: 

PO «2 - PPO COS 61 [62 - 4''] + /2 = 0, 

PO 02 + PPo COS 61 [62 - cl^] + /2 + PPo 02 = 0, 

witfi 

f2{x) = ?72 cos^ 9 d^h-i + (sin 6* - 9^/i_i)9^cf . 

Concentration: 

2 /2(:r) 



02 (x) = 



Po{x) 2po{x) +p' 
02 (x) 



2cos^' 

1 \ . . 

x). 



C2ix,z) = 4\^)+{z-j^) a2 
Zero normal velocity. 

order 3 

Diffusion equation: 

Po (^zzC3 - 2pl ds = 0. 

with 

. s _ 2(a^/i_i - sin 6*) 9^0,2 + 0,2 d^^h_i - d^^cl^ 

Solution: 

C3(x, z) = ci3(a;)-s^ + 03(0;) 2; + b3{x). 



Boundary conditions: 



pI 03 - PPo cos 9 [bs - + /3 = 0, 

pI 03 + PPo cos 9 [63 - C3''] + /3 + 5-3 + PPo 03 = 0, 
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with 



h[x) = d^a2 - o^ho d^c^* 

2 cost^ 

2{dxh_i- sm.9] d^h^ - p^pi 



+ 



02 



+ (c^ - — ) dxh_i + dnchQ r]2 cos^ 6 
+ {dxh^i - sin 9^ {p- d^c^ - d^c. 



2 /' 



and 



Concentration: 



Po(2j9o + p) , - sin6l 
S'aW = —7, 03 —n o^a2. 



cos Q 



cos 6 



(79) 



(80) 



03 (X) 



± 93 + 2/3 
Po 2po + p 



p/3 - Pogs 



ppo COS (2po + p) 
(i3(x)z^ + a3(x)z + 63(2;) 



(81) 



The normal velocity is nonzero for the first time at this order. Its expression is derived by 
using Eq. (jl^ together with the scaling relations of section IIVC 1\ 



a 



2 COS^ + PPO r^ eq 

'Ox / ^ X OxC^ 



Pqcos'^0 [ po(p + 2po) 

2 /) sin 6* - dxh_i 

- 2r/2 cos — 

Po(p + 2po) 



(82) 



D. Amplitude equation 



We finally obtain the following amplitude equation using Eqs. ()49I82|) : 



dtH 



a 



cos^ 



2cOs2 6' + PP0q ^eq 



_ po{p + 2po) 



-dxC^^ 



2/1 sine - dxH 
2?72COS : z — rO^.H 



where 



and 



Po(p + 2po) 
Poix) = ^{1 - sm0 dxH)^ + {cos0 dxH)^, 



Ci (x) = —acos0 — 5— 



l3 + apl[ 2pl - cos^ 



(83) 



(84) 



(85) 
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FIG. 6: Numerical simulation of Eq. (|83j) . Time evolution of a single step for p = 0.001 and 
9 = 0.8. The step is systematically shifted in time (given by the lower axis). The electrical field is 
applied in the positive x direction. 

Here H{x) = h-i{x) and the time is rescaled such as eH t. This amplitude equation is 
the central result of our study. As expected, this equation ensures mass conservation since 
its right hand side is a derivative of a mass current. 

V. NUMERICAL SIMULATIONS AND DISCUSSION 

The time evolution of vicinal surfaces is obtained by integrating numerically Eq. ()83p . 
While the simulations are performed in the oblique frame (a;,z), the system is represented 
in the laboratory orthogonal frame {x,y). Solving this stiff partial differential equation ne- 
cessitates the use of an adaptive time step. A single step with periodic boundary conditions 
is simulated in practice. The whole vicinal surface is obtained by reproducing this step 
periodically along the z direction. The elastic interactions included in our model are not 
only justified from a purely physical point of view but are also a necessary ingredient in real- 
istic numerical simulations. Indeed, test simulations performed without elastic interactions 
systematically resulted in step crossings at late times. 

We first compare the dynamics of one step in two physical regimes defined by the values 
of the nondimensional number p = For p > 1, the system dynamics is diffusion-limited. 
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FIG. 7: Same as in Fig. for p = 20 and a larger system width. 

while it is attachment-limited for p < 1. All the parameters {a, f3, 772, a) entering Eq. ()83|) 
are set to unity here, and Figs. © and (|7j) show the time evolution of a single step for 
p = 0.001 and p = 20, respectively. At short times, the steps are rather similar in shape 
for both values of p. Calculating the wave length emerging at short times, we find that it 
increases with p as predicted by the linear stability analysis. Alternatively, the growth rate 
r is found to decrease with p. At late times, after coarsening has set in, the step shapes differ 
strongly: a single-valued function is found in the laboratory frame for p = 0.001, while long 
overhangs are visible for p = 20. In both cases, the electrical field triggers local facetting 
of the steps which look like asymmetrical saw-teeth. Ultimately, the meander amplitude 
saturates to a finite value in a finite size system. 

The time evolution of two vicinal surfaces is displayed in Figs.® and Dark regions 
correspond to a high step density in which the electrical field is essentially oriented in the 
step-down direction, while it is mainly oriented in the step-up direction in the low step 
density regions. This result is consistent with the well-known step bunching observed for 
Si (111) when the heating current is applied perpendicular to the steps, in the step-down 
directioniii. 

According to the nature of the material, the surface orientation and the temperature 
range, physical parameters such as the diffusion coefficient may vary a lot. In addition. 
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(a) (b) 




(c) (d) 

FIG. 8: Top view of a vicinal surface computed at different times for the same parameters as in 
Fig. ©: a) t = 230, b) t = 1200, c) t = 8300, d) t = 1.75 x 10^. The step down direction is 
rigthwards while the electrical field direction is downwards. 

they are not always known with a great accuracy. For example, for a Si(lll) surface, four 
acceptable sets of physical parameters are given in table I of ref-^^, of which set B seems 
particularly consistent with the experimental observations. For this particular set of physical 
parameters, Eq. (fT^ gives d = Lq/ p = 5x 10^ ''m. A miscut angle of one degree, then results 
in p ~ 0.03, thus attachment/detachment-limited dynamics. Note that with the parameter 
sets A, C, or D, and/or a different miscut angle, p may vary in wide range, both below and 
above one. Our model is valid in both cases and it predicts rather different step shapes at 
long times, as just discussed. Experimental observation of vicinal surfaces under an electrical 
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FIG. 9: Top view of a vicinal surface computed at different times for the same parameters as in 
Fig. Q: a) t = 6.4 x 10^, b) t = 1.6 x 10^ c) t = 8 x 10^ d) t = 3 x 10^ The step down direction 
is rigthwards while the electrical field direction is downwards. 

field parallel to the initial steps could possibly give an indication on the magnitude of the 
nondimensional number p which governs the system dynamics. 

VI. CONCLUSION AND PERSPECTIVES 

In summary, we have studied the meandering instability induced by a constant electrical 
field initially parallel to a train of straight steps. The time evolution of the meanders is 
described by a nonlinear amplitude equation which we have derived through an asymptotic 
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expansion. Numerical simulations have been performed both in the attachment/detachment- 
limited (p <^ 1) and the diffusion- limited (p ^ 1) regimes. At large times, overhangs are 
observed in the latter case only. 

It is very instructive to compare our results with an experimental study of step me- 
andering on Si (111) vicinal surfaces, in which the orientation of the electrical field E is 
taken different from the step-down direction^''. When E is set parallel to the steps, as in the 
present study, a similar step meandering effect is observed but the steps bend in the opposite 
direction as compared to our model. This apparent contradiction is in fact not unexpected 
because the experiments are performed at T = 1100°C. Indeed, in this intermediate range 
of temperature (1000°C — 1180°C), the steps have been argued to become transparent to 
the diffusing adatoms^*^. The underlying physics is thus expected to differ from the one 
introduced in our model (impermeable steps) and an opposite direction of bending is not 
contradictory. In the light of this discussion, new experiments performed at temperatures 
slightly higher than T = 1180°C or slightly lower than T = 1000°C would be desirable to 
test our model. 

In the present model, consecutive steps are assumed identical up to a given phase-shift. 
Removing this phase constraint would allow a realistic description of experiments on a large 
scale. However, this can hardly be envisaged on the basis of the present method and a 
quite different point of view should be considered, such as a continuous limit approach. 
In addition, it would be helpful to include the step transparency in order to compare the 
resulting model to the experiments in the intermediate range of temperatures. 
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